from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from numpy import arange


def f(T, H, SC, k_1, k_3, k_4, k_5, k_6, k_H, S_C, S_S, eta0, n, alpha):
    # 变量
    solid_content = SC / 100  # 固含量
    T = T + 273.15  # 温度，单位为开尔文

    # 常量
    C_0 = 24e-3  # 比例系数，用于确定DMF质量与时间的关系
    m_D_0 = 24e-3  # 初始时刻DMF的质量，单位为kg
    m_S = 6e-3  # 环丁砜总的质量，单位为kg。无论是否溶解，环丁砜不会从体系中消失，总的质量是不变的
    m_C = (m_D_0 + m_S) * solid_content  # 初始时刻醋酸纤维素的质量，单位为kg，与初始固含量有关
    ro_D = 0.948e3  # DMF的密度，单位为kg/m^3
    ro_S = 1.261e3  # 环丁砜的密度，单位为kg/m^3
    ro_C = 1.3e3  # 醋酸纤维素的密度，单位为kg/m^3
    V = m_D_0 / ro_D + m_S / ro_S + m_C / ro_C  # 总体积，单位为m^3，认为不变
    A = 6.09451
    B = 2725.96
    C = 28.209  # 由温度计算饱和汽压时的三个常数
    A_0 = -k_1 * (10 ** (A - B / (T + C))) * 133.322 * (1 - k_H * H / 100) / ro_D / V
    n_correction_factor = n  # 斯托克斯 - 爱因斯坦方程中的修正系数，对于大分子乃至宏观颗粒是6，对于小分子则体积越小，该值也会越小
    k = 1.380649e-23  # 玻尔兹曼常数，单位为J/K
    r = 0.3413e-9  # 环丁砜的分子半径
    N_A = 6.02214076e23  # 阿伏伽德罗常量，单位为mol^-1

    # 函数
    # DMF在蒸发的过程中，其质量随时间变化的函数
    def m_D(t):
        return C_0 * np.exp(A_0 * t)

    # 未能溶解而析出成液滴的环丁砜的质量随时间变化的函数
    def m_S_out(t):
        return max(0, m_S - m_D(t) * S_S)

    # 未能溶解而析出成固体的醋酸纤维素的质量随时间变化的函数
    def m_C_out(t):
        return max(0, m_C - m_D(t) * S_C)

    # 析出的醋酸纤维素的体积占比随时间变化的函数
    def phi_C_out(t):
        return m_C_out(t) / ro_C / V

    # 不同温度下混合溶液总的粘度随时间变化的函数
    def eta(temp, t):
        return eta0 * (1 + 2.5 * phi_C_out(t))

    # 扩散系数。假定一个小液滴内仅有一个分子，且分子为球形
    def D(temp, t):
        return k * temp / (n_correction_factor * np.pi * r * eta(temp, t))

    # 液滴运动的平均速度
    def v(temp, t):
        return k_3 * (D(temp, t) ** 0.5)

    # 析出的环丁砜的分子数密度随时间变化的函数
    def n_molecular_number_density(t):
        return m_S_out(t) / (ro_S * 4 / 3 * np.pi * r ** 3) / V

    # 小液滴间的碰撞频率
    def Z(temp, t):
        return (2 ** 0.5) * n_molecular_number_density(t) * np.pi * (r ** 2) * v(temp, t)

    tf = np.log(alpha) / A_0  # 对于每种情况，在反应了一定时间后停止

    dy = lambda m_s, t: k_4 * Z(T, t) + k_5 * v(T, t) * (m_S_out(t) - m_s) / V * m_s
    t = arange(1, tf, 0.01)
    sol = odeint(dy, 0, t)

    if len(sol.T[0]) == 0:
        re = 10000
    else:
        re = k_6 * sol.T[0][-1]
    return re


bound = np.array([
    [30, 50],
    [50, 90],
    [6, 10],
])


# 在SA函数外部初始化记录容器
history = {
    'temperature': [],
    'best_loss': [],
    'current_loss': [],
    'datas': []
}

data_name = ["温度", "湿度", "固含量"]


def plot_SA_history():
    plt.figure(figsize=(14, 10))
    plt.rc('font', size=16)
    plt.rc('font', family='SimHei')

    # 子图1：损失函数变化
    plt.subplot(2, 2, 1)
    plt.semilogy(history['best_loss'], 'r-', label='Best Loss')
    plt.semilogy(history['current_loss'], 'b--', alpha=0.5, label='Current Loss')
    plt.xlabel('Iteration')
    plt.ylabel('Loss')
    plt.title('Loss变化曲线')
    plt.legend()
    plt.grid(True, which="both", ls="--")

    # 子图2：温度衰减曲线
    plt.subplot(2, 2, 2)
    plt.plot(history['temperature'], 'g-')
    plt.xlabel('Iteration')
    plt.ylabel('Temperature')
    plt.title('温度下降曲线')
    plt.grid(True, ls="--")

    # 子图3：k参数演化
    plt.subplot(2, 2, 3)
    k_array = np.array(history['datas'])
    for i in range(3):
        plt.plot(k_array[:, i], label=f'{data_name[i]}')
    plt.xlabel('Iteration')
    plt.ylabel('温度、湿度、固含量')
    plt.title('制备条件变化曲线')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, ls="--")

    plt.tight_layout()
    plt.show()


def SA(t0, tf, alpha, iter):
    t = t0
    data0 = np.array([40, 70, 8])
    datac = data0
    datab = data0
    lc = f(40, 50, 6, 0.0010660234370888787, 0.34052550181423463, 0.11440264521150946, 0.03446581365646933, 0.30328518806514704, 0.32248515839275166, 0.9118799744432465, 0.6509464523758399, 0.700032233700628, 2.495457103926549, 0.2574903993186905)
    lb = f(40, 50, 6, 0.0010660234370888787, 0.34052550181423463, 0.11440264521150946, 0.03446581365646933, 0.30328518806514704, 0.32248515839275166, 0.9118799744432465, 0.6509464523758399, 0.700032233700628, 2.495457103926549, 0.2574903993186905)
    for i in range(iter):
        datan = datac + np.random.normal(0, 0.5, size=3)
        datan = np.clip(datan, bound[:, 0], bound[:, 1])
        ln = f(datan[0], datan[1], datan[2], 0.0010660234370888787, 0.34052550181423463, 0.11440264521150946, 0.03446581365646933, 0.30328518806514704, 0.32248515839275166, 0.9118799744432465, 0.6509464523758399, 0.700032233700628, 2.495457103926549, 0.2574903993186905)

        if ln > lc or np.random.rand() < np.exp((ln - lc) / t):
            datac = datan
            lc = ln
            if lc > lb:
                datab = datac
                lb = lc
        t *= alpha
        if t < tf:
            break
        print(f"iter{i}, lb{lb},datac{datac},lc{lc},dets_loss{ln-lc}")
        history['temperature'].append(t)
        history['best_loss'].append(lb)
        history['current_loss'].append(lc)
        history['datas'].append(datac.copy())
    print(f"最优孔面积占比：{lb}")
    return datab


re = SA(300, 0.001, 0.99, 1000)
plot_SA_history()
print(re)
